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The  multiple-relaxation-time  (MRT)  lattice  Boltzmann  method  (LBM)  with  multi-reflection  solid  bound¬ 
ary  conditions  is  used  to  study  anisotropic  permeabilities  of  a  carbon  paper  gas  diffusion  layer  (GDL)  in  a 
fuel  cell.  The  carbon  paper  is  reconstructed  using  the  stochastic  method,  in  which  various  porosities  and 
microstructures  are  achieved  to  simulate  different  samples.  The  simulated  permeability  and  tortuosity 
show  anisotropic  characteristics  of  the  reconstructed  carbon  papers  with  in-plane  permeability  higher 
than  through-plane,  and  in-plane  tortuosity  lower  than  through-plane.  The  calculated  permeabilities  are 
in  good  agreement  with  existing  measurements.  The  relationship  between  the  permeability  and  the  poros¬ 
ity  is  fitted  with  empirical  relations  and  some  fitting  constants  are  determined.  Furthermore,  the  obtained 
relationship  of  tortuosity  and  porosity  is  used  in  a  fractal  model  for  permeabilities.  The  results  indicate 
that  the  fractal  model  and  the  Kozeny-Carman  equation  provide  similar  predictions  on  the  through-plane 
permeability  of  the  carbon  paper  GDL. 

©  2008  Published  by  Elsevier  B.V. 


1.  Introduction 

Due  to  their  high  efficiency,  low  pollution  and  low  noise,  pro¬ 
ton  exchange  membrane  fuel  cells  (PEMFCs)  have  attracted  much 
attention  for  over  a  decade  as  a  promising  candidate  for  power 
sources  in  automotive  and  other  portable  electronic  devices.  A  typ¬ 
ical  PEMFC  consists  of  a  polymer  membrane  sandwiched  between 
two  electrodes  (anode  and  cathode).  Each  electrode  can  be  divided 
into  three  regions:  the  bipolar  plate,  the  gas  diffusion  layer  (GDL) 
and  the  catalyst  layer.  In  PEMFCs,  fuel  and  oxidant  are  supplied 
to  the  channels  on  the  bipolar  plate,  and  transfer  through  the  gas 
diffusion  layer  to  catalyst  activity  sites  where  electrochemical  reac¬ 
tions  occur.  The  GDL  plays  an  important  role  for  providing  structural 
support,  permeating  reactant  gas,  removing  product  water  and 
conducting  electrons.  The  permeation  of  gas  and  removal  of  liq¬ 
uid  water  from  GDL  are  critical  to  the  performance  of  the  PEMFCs, 
because  higher  reactant  mass  transfer  rates  would  achieve  higher 
current  densities  while  the  accumulation  of  liquid  water  in  GDL 
would  block  reactant  transfer,  thus  lowering  the  performance  [1]. 
Therefore,  the  permeability  of  the  GDL  is  a  key  parameter  on  the 
performance  of  the  PEMFCs. 
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Generally,  the  GDL  is  a  porous  medium  made  of  carbon  paper  or 
carbon  cloth.  These  materials  have  microscopically  complex  struc¬ 
ture  with  random  distribution  of  pore  sizes,  ranging  from  a  few 
microns  to  tens  of  microns.  In  practice,  in  order  to  decrease  the 
resistance  for  the  electrons  transfer  from  the  catalyst  layer  to  the 
bipolar  plane,  the  GDL  is  slightly  compressed  by  external  force, 
by  which  the  permeability  is  reduced  [2].  Experimental  studies 
confirm  that  there  exists  an  optimal  compression  for  the  GDL, 
attributing  to  the  trade-off  between  the  improved  contact  resis¬ 
tance  and  the  reduced  GDL  permeability  [3].  Additionally,  it  is 
believed  that  treating  the  GDL  with  a  hydrophobic  polymer  such  as 
PTFE  would  improve  water  removal  and  thus  facilitate  gas  diffusion 
in  the  GDL  However,  studies  also  show  that  electrical  conductiv¬ 
ity  and  permeability  of  the  GDL  are  reduced  if  excessive  PTFE  is 
added  [4],  due  to  the  insulating  property  of  PTFE  and  the  filling  of 
the  GDL  pores.  Thus,  finding  optimal  values  of  the  GDL  parameters 
such  as  porosity,  permeability  and  wettability  are  important  issues 
in  PEMFCs  research  [5-7]. 

Recently,  a  great  deal  of  attention  is  given  to  the  study  of 
anisotropic  characteristics  of  the  GDL,  i.e.,  permeabilities  are 
different  in  the  in-plane  and  in  the  through-plane  directions.  Mea¬ 
surements  of  through-plane  permeability  were  performed  [7,8], 
and  the  correlation  between  the  through-plane  permeability  and 
the  limiting  current  density  was  studied  [9].  In  addition,  several 
studies  [10-12]  indicated  that  in-plane  permeability  is  more  rele¬ 
vant  to  PEMFCs  performance,  particularly  for  the  serpentine  flow 
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channel  in  which  the  cross-flow  occurs  due  to  pressure  gradient 
between  the  neighbor  channels.  Feser  et  al.  [13]  and  Gostick  et  al. 
[14]  measured  in-plane  permeability  of  carbon  paper  and  carbon 
cloth  GDL  as  a  function  of  compression  ratio,  and  found  that  the 
permeability  of  the  GDL  decreased  with  the  increase  of  the  com¬ 
pression  force. 

Early  work  on  simulation  of  fluid  flow  in  the  GDL  used  macro¬ 
scopic  models  based  on  the  volumetric  averaging  of  conservation 
equations  in  an  elementary  volume  assuming  a  homogeneous 
porous  medium  with  a  constant  permeability,  which  were  solved 
numerically  by  conventional  CFD  methods  [15].  Recently,  some 
work  on  the  application  of  the  lattice  Boltzmann  method  (LBM) 
to  study  transport  in  GDL  of  PEMFCs  have  been  carried  out  [16-18]. 
Using  single-phase  LBM  approach,  Wang  and  Afsharpoya  [16]  stud¬ 
ied  the  2D  fluid  flow  though  a  section  of  serpentine  channel  and  a 
channel  filled  or  partially  filled  with  porous  medium  with  appli¬ 
cation  to  PEMFCs  gas  flow  channel.  Joshi  et  al.  [17]  used  the  LBM 
to  model  the  mass  transport  of  H2  and  the  produced  H20  (vapor) 
in  the  presence  of  N2  in  the  2D  porous  anode  structure  of  a  solid 
oxide  fuel  cell  (SOFC).  Park  and  Matsubara  [18]  studied  the  gas 
flow  through  a  fiber  tow  in  the  woven  carbon  cloth  GDL  by  using 
a  3D  LBM  approach  in  solving  the  Stokes  and  Brinkman  equations 
numerically  in  the  void  and  in  the  porous  tow,  respectively.  The  LBM 
used  in  [18]  is  based  on  the  Bhatnagar-Gross-Krook  (BGK)  model 
with  a  standard  bounce-back  boundary  scheme  for  solid  bound¬ 
ary  condition.  The  BGK  model,  which  has  a  single-relaxation-time 
(SRT)  collision  operator,  suffers  from  numerical  instability  and  vis¬ 
cosity  dependence  of  flow  when  solid  boundary  is  present  [19]. 
The  viscosity  dependence  problem  becomes  severe  for  simulating 
flow  through  porous  media.  To  overcome  these  deficiencies,  the 
multiple-relaxation-time  (MRT)  model  has  recently  been  devel¬ 
oped  [19-21],  in  which  different  relaxation  times  are  used  for 
different  kinetic  modes. 

In  this  paper,  the  3D  micro-scale  flow  of  liquid  transport  in 
the  nonwoven  carbon  paper  GDL  is  simulated  using  the  MRT/LBM 
approach  [19,20].  In-plane  and  through-plane  permeabilities  of  the 
anisotropic  GDL  are  simulated  for  several  samples  that  are  recon¬ 
structed  using  the  stochastic  method  [22,23].  The  samples  with 
different  porosities  and  microstructures  are  achieved  to  simulate 
different  GDL  compression  and  PTFE  contents  and  the  calculated 
permeabilities  are  compared  with  existing  measurements  [  13 ].  The 
numerically  obtained  relationship  between  the  permeability  and 


the  porosity  is  compared  with  the  previously  developed  correla¬ 
tions  and  with  fractal  models  [24,25]. 


2.  Description  of  the  simulation  model 


2.1.  Multiple-relaxation-time  LBM  model 


The  MRT/LBM  model  transforms  the  distribution  function  in 
the  velocity  space  of  the  SRT/LBM  model  to  the  moment  space 
through  a  transformation  matrix.  Since  moments  of  the  distribu¬ 
tion  function  directly  represent  physical  quantities,  the  moment 
representation  offers  a  convenient  way  to  perform  the  relaxation 
collision  processes  with  different  relaxation  times  according  to  dif¬ 
ferent  time  scales  of  various  physical  processes.  In  the  SRT/LBM 
model,  the  evolution  of  the  distribution  function/a(x,  t)  is  described 
by  the  following  equation: 

fa{r+ea8t ,  t+<$t)-/„(r,  t)=S  [ faeq(r+ea8t ,  t+8t)-fa( r,  f )]  +Fa  (1 ) 

where  S  is  the  relaxation  matrix  [21],  which  is  a  diagonal  matrix 
with  same  diagonal  elements  in  SRT  model;  ea  is  lattice  velocity  in 
the  ath  direction  and  the  corresponding  lattice  velocity  directions 
are  shown  in  Fig.  1  .f^q  is  the  equilibrium  distribution  of fa  given  as 
a  function  of  density  p  and  velocity  u: 


faq  =  Wap 


(eg  •  u)2 
2cs4 


u  u 


(2) 


For  the  three-dimensions  19-velocity  (3DQ19)  model,  the  lattice 
velocity  ea  and  the  weight  coefficients  wa  are  given  as  follows: 


^cy  — 


Wa  = 


(0,  0,  0),  <2=0; 

(±1,0,  0)c,  (0,±l,0)c,  (0,0±l)c,  <2= 1,2, 

(±1,  ±1,  0)c,  (±1, 0,  ±l)c,  (0,  ±1,  ±l)c,  <*=7,8, 
1/3,  a  =  0; 

1/18,  a  =  l,2,  ...,6; 

1/36,  <2  =  7,8, ...,  18. 


6; 

18. 


(3) 


Fa  in  Eq.  (1 )  represents  the  external  force,  which  is  given  by  [26] 


(4) 


el2 


Fig.  1.  Lattice  velocity  directions  of  the  D3Q19  LBM. 


106 


L.  Hao,  P.  Cheng  /  Journal  of  Power  Sources  186  (2009)  104-114 


where  F  is  the  external  force  and  cs  is  the  speed  of  sound  which 
is  given  by  c/V 3,  where  c  =  8x/8t  is  the  lattice  speed  and  8x  is  the 
lattice  distance. 

By  multiplying  a  transformation  matrix  Q  in  Eq.  (1),  the  evolu¬ 
tion  function  can  be  expressed  in  the  moment  space  as 

m(r+e<5t,  t+8t)-m{r,  t)= S  [meq(r+e^t,  t+8t)- m(r,  t)]  +Q  F 

(5a) 

where 

m  =  Q  f  (5b) 


Employing  the  Chapman-Enskog  multiscale  analysis  for  D3Q19 
model,  Eqs.  (1)  and  (5)  can  lead  to  the  Navier-Stokes  equations 
under  the  low  Mach  number  limitation: 


^  +  V  ■  (pu)  =  0 

(10a) 

9(9fU)  +  V  ' (/0UU)  =  V(Cs2/0)  +  V  '  ^yVu)  +  F 

(10b) 

where  the  pressure  is  given  by  p  =  cjp. 

2.2.  Boundary  conditions 

and 

m  e«  =  (l.fq  (5c) 

are  the  vectors  of  the  moments  distribution  function  and  its 
corresponding  equilibrium  moments  vector,  respectively.  The 
transformation  matrix  Q  is  constructed  such  that  the  relaxation 
matrix  S  in  moment  space  can  be  reduced  to  the  diagonal  form  and 
is  given  in  [20]  for  D3Q19  model.  As  a  result,  the  diagonal  relaxation 
matrix  S  in  Eq.  (5a)  is 

S  =  diag(So,  Si ,  S2,  S3,  S4,  S5,  S6,  S7,  S8,  Sg,  Sio,  Sn  ,  Si2,  Si3,  Si4, 

Sl5>Si6, S17,  Sis)  (6) 

where  s0  =s3  =Ss  =S7  =  0,  Si  =s2  =s9_i5  =  1/r,  and  s4  =  s6  =  s8  = 
$16-18  =  8((2  -  1  /r)/(8  -  1  /r))  are  propositional  values  in  [19]  and 
r  is  related  to  the  fluid  viscosity  v  by 


Additionally,  the  equilibrium  moments  meq  are  [19]: 
nC  =  p  (8a) 


Since  the  methods  developed  in  the  SRT  model  for  boundary 
conditions  are  still  applicable  in  the  MRT  approach,  the  moment 
distribution  function  m  in  the  MRT  model  is  transformed  to  the 
distribution  function  f  in  the  velocity  space  by 

f=QT1m  (11) 

when  dealing  with  the  boundary  conditions. 

2.2.1.  Inlet/ outlet  boundary 

The  velocity  and  pressure  boundary  conditions  at  the  inlet  and 
outlet  are  dealt  with  using  the  bounce-back  of  the  non-equilibrium 
distribution  rule  [27].  For  the  D3Q19  model,  for  example,  if  the  inlet 
boundary  face  is  perpendicular  to  the  x-direction  with  the  lattice 
velocity  e\ ,  e7,  ei0,  en  and  ei4  (as  shown  in  Fig.  1 )  pointing  into  the 
calculation  region,  the  distribution  functions  f\ ,  /7,  /i0,  fu  and  /i4 
in  these  directions  are  unknown  after  streaming.  These  unknown 
distribution  functions  are  determined  by  the  rest  of  the  known 
distribution  functions  in  order  to  satisfy  the  specified  boundary 
conditions. 

According  to  Eqs.  (9a)  and  (9b),  the  unknown  distribution  func¬ 
tion//^  1,  7, 10, 11, 14)  at  the  inlet  boundary  satisfy  the  equations 
as  follow: 


m 


eq 

1 


-llp  + 


19  — 
Po 


meq 

m3 

—  jx, 

meq 

m5 

=  jy , 

3 

vj  to 
-Q 

II 

nT‘ 

me4q 

2. 

3  Jx  1 

2  .  eq 

=  -3 Jy’  ms 

2. 

-  -3  Jz 

m9 

 3 &  -j 

Po 

J 

meq  3j2  -  j  j 

10  “  2p0 

mn 

_  jy  -iz 
Po 

’ 

mi2 

Jy-Jz 

2po 

< 

_jxjy 

Po  ’ 

meq 

mu  = 

4k,  = 

Po  15 

jxjz 

Po 

mi6 

=  mi7  = 

< 

=  0 

(8b) 

(8c) 

(8d) 

(8e) 

(8f) 

(8g) 

(8h) 

(8i) 


f\  +fl  +/l0  +/ll  +/l4  =  Pin  ~  ifo  +/2  +/3  +/4  +/s  +  fo  +/8  +  k 

+/12 +/13 +/15 +/16 +/17 +/18)  (12a) 

/l  T  fl  +/l0  +/ll  +/l4  =  Pin  —  +  C/4  +/8  T /9  +/12  +/13)  (12b) 

where  the  inlet  density  pin  can  be  obtained  by  equating  the  right- 
hand  side  of  Eq.  (12a)  to  the  right-hand  side  of  Eq.  (12b)  to  give: 


E 

Pin  ~  T~^/-c 


(13a) 


where 

E  =  ifo  +/2  +h  +/s  +fe  +/15  +/16  +/17  +/18  +  2  Cf4  +/8  +/9 

+/12+/13)  (13b) 


To  close  the  system  of  equations,  the  bounce-back  rule  for  the 
non-equilibrium  part  of  the  distribution  functions  fi(i  =  1,  7, 10,  11, 
14)  is  applied 


where  po  is  the  mean  density  of  the  fluid  which  is  used  to  reduce 
compressibility  effects  in  the  model  [19,20]. 

The  macroscopic  local  density  p(x,  t)  and  the  momentum  j(x,  t) 
are  obtained  as 
18 

p = yya  (9a) 

a=0 

+  yF  (9b) 

a=0 


j  =  p0U  =  Yjfc 


+  (14) 

with  ffq  calculated  from  pin  given  by  Eq.  (13)  and  with  f  denotes 
the  opposite  lattice  velocity  of  i,  that  is  ez  =  -ev.  Eqs.  (12)  and  (14) 
are  sufficient  to  determine  all  the  unknown  //.  However,  in  order 
to  keep  the  correct  y-  and  z-direction  momenta,  these  distribution 
functions  are  modified  as  follows  [27]: 

^  =  fi-Uyei^+jzeiz)  ^  i  =  i7W,n,U  (15) 

where  jy  =  J2li(/«e«v  andJz  =  Eaii /“e“z- 
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Similarly,  when  the  boundary  pressure  (density)  is  specified,  the 
x-direction  velocity  can  also  be  determined  by  Eq.  (12),  i.e., 


ux  =  c 


(16) 


than  bounce-back  rule,  especially  for  porous  media  flow  [19,20].  In 
addition,  the  MR  boundary  condition  utilizes  some  moment  space 
parameters,  which  can  be  directly  obtained  in  the  MRT/LBM  model. 

After  the  streaming  process,  the  distribution  function  out  of  the 
solid  is  constructed  in  MR  boundary  condition  as  follows: 


Under  the  assumption  that  velocities  in  y-  and  z-directions  are 
zero  on  the  boundary,  the  unknown  distribution  functions /j(i  =  1, 
7, 10, 11, 14)  can  be  obtained  using  Eqs.  (14)  and  (15)  as  the  same  as 
velocity  boundary. 


fr{x,  t )  =  kifi(x  +  eiSt,  t )  +  krfifc  t )  +  k3fi(x  -  efit,  t )  +  krf^x,  t) 
+ksfi'(x  -  e,St,  t)  +  k6Ny(x,  t)/v  (17) 


2.2.2.  Solid  boundary 

The  no-slip  condition  on  the  solid  boundary  can  be  approxi¬ 
mated  using  the  standard  bounce-back  boundary  conditions,  in 
which  the  liquid  particles  reflect  their  directions  when  colliding 
with  the  solid  wall.  However,  it  has  been  demonstrated  that  the 
multi-reflection  (MR)  boundary,  which  was  developed  originally 
for  solving  moving  solid  boundary  problems,  is  much  more  accurate 


with  the  coefficients  kj  (j  =  1  -6)  being  functions  of  the  solid  bound¬ 
ary  location  given  by  [19] 


k\  =  1 ,  k2  =  -k4  = 


1-2  q-  2  q2 
(1  +<7)2 


k3  =  -k5  = 


kfi  = 


1 


4(1  +qf 


(i  +  <?r 


Fig.  2.  Micrographs  of  the  carbon  paper  GDL.  (a)  SEM  image  of  the  Toray  060  GDL.  (b)  3D  geometer  of  the  reconstructed  carbon  paper  GDL.  (c)  3D  geometer  of  the  reconstructed 
carbon  paper  GDL  with  15%  PTFE  content. 
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where  q  is  defined  as  the  ratio  of  the  distance  between  the  first 
fluid  node  near  the  wall  and  the  solid  wall  to  the  lattice  distance. 
In  Eq.  (17)  the  quantity  N  is  the  non-equilibrium  term  relating  to 
moments  distribution  functions  m  and  can  be  calculated  as 

N  =  M-1  •  S  •  (tc?  —  t)  (18a) 

where  t  =  (0,  0,  0, 0,  014,  0,  mg,  0,  mg,  0, 0,  m^g,  m^,  mis)T  is  a 
vector,  having  partial  elements  of  the  moments  m,  and 

teq  =  (0,  0,  0,  0,  meq,  0,  meq,  0,  me8q,  0, 0,  meq6,  m eqn,  m^)T  (18b) 

The  main  steps,  when  applying  the  MRT/LBM  model  described 
above,  can  be  summarized  as  follows: 

1.  Calculating  the  equilibrium  moments  meq  from  Eq.  (8); 

2.  Evolution  of  the  moments  distribution  function  according  to  Eq. 

(5); 

3.  Transforming  the  moment  distribution  function  m  to  the  veloc¬ 
ity  distribution  function  f  by  Eq.  (11 ); 

4.  Implementing  the  streaming  process  using  the  velocity  distribu¬ 
tion  function  f  and  imposing  the  boundary  condition  from  Eqs. 
(12)— (18); 

5.  Updating  the  local  density  p(x,  t)  and  the  local  velocity  u(x,  t) 
according  to  Eq.  (9). 


original  samples  and  does  not  change  the  fibers  shape.  Thus,  the 
reconstruction  process  is  the  same  as  mentioned  above,  only  with 
different  porosities  and  thicknesses  according  to  different  compres¬ 
sion  ratios.  The  relationship  between  the  compression  ratio  r  and 
the  porosity  0  can  be  expressed  as  [13] 

0  =  1-  t~t  (19) 

where  0O  is  the  uncompressed  GDL  porosity,  and  r=(h0  -  h)/ho  is 
defined  as  the  change  of  the  thickness  h0-h  after  compression 
divided  by  the  uncompressed  thickness  h0.  Thus,  r=0  means  no 
compression. 

If  the  PTFE  is  added  on  the  carbon  paper  GDL,  one  more  step  is 
added  after  the  fiber  GDL  generation:  the  pore  volume  close  to  the 
fiber  in  the  sample  is  randomly  filled  by  the  PTFE.  This  process  is 
repeated  until  the  prescribed  PTFE  content  is  achieved.  Thus,  the 
pores  surrounding  by  more  fibers  will  have  more  probability  to  fill 
with  PTFE.  Fig.  2(c)  shows  the  geometry  of  the  reconstructed  GDL 
with  15%  PTFE  content  and  the  same  size  of  150  x  150  x  127  lat¬ 
tice  grid.  The  final  porosity  0  can  be  expressed  approximately  as  a 
function  of  PTFE  content  w  to  give: 

,  ,  w(l-0o) 

4>  =  4>o-a— j— —  (20) 


2.3.  Reconstruction  of  carbon  paper  GDL 


where  0O  is  the  original  porosity  before  PTFE  is  added,  and  a  =  0.9 
is  the  density  ratio  of  the  carbon  fiber  and  the  PTFE  [7]. 


The  reconstruction  of  3D  porous  media  can  be  obtained  by  two 
methods:  the  3D  imaging  combination  and  the  virtual  stochas¬ 
tic  generation.  The  former  employs  the  3D  images  of  the  porous 
medium  scanned  by  X-ray  or  scanning  laser  microscopy  and  then 
integrates  them  to  achieve  the  geometry  construction.  The  later 
reconstructs  the  porous  medium  microstructure  based  on  the  sta¬ 
tistical  information  of  the  porous  medium.  The  low  cost  and  easy 
implementation  of  geometry  generation  make  the  stochastic  gener¬ 
ation  method  a  better  choice  than  the  imaging  combination  method 
[22].  In  the  present  study,  the  method  proposed  by  Schladitz  et 
al.  [23]  based  on  the  stochastic  generation  method  is  adopted  to 
reconstruct  the  carbon  paper  GDL. 

Fig.  2(a)  presents  the  SEM  image  of  the  Toray  060  carbon  paper. 
In  general,  an  actual  carbon  paper  GDL  consists  of  carbon  fibers, 
randomly  oriented  in  a  plane,  leading  to  different  in-plane  and 
through-plane  properties.  Since  the  reconstruction  of  the  real  GDL 
is  difficult,  some  assumptions  are  made  for  simplification:  (i)  car¬ 
bon  fibers  in  GDL  are  cylinders  having  a  fixed  and  uniform  diameter; 
(ii)  fibers  are  straight  and  infinitely  long;  (iii)  fibers  are  allowed 
to  overlap.  In  addition,  according  to  fabrication  process  of  carbon 
paper,  the  fibers  in  the  material  plane  are  stochastic  arranged,  thus 
this  plane  can  be  considered  to  be  isotropic,  and  fibers  oriented 
to  the  direction  of  the  GDL  thickness  can  be  neglected.  Giving  the 
GDL  size,  porosity  and  fiber  diameter  D,  the  reconstruction  of  car¬ 
bon  paper  is  achieved  by  placing  the  fibers  in  a  plane  with  arbitrary 
positions  and  orientations  one  layer  by  one  layer.  It  should  be  noted 
that  the  porosity  of  the  reconstructed  GDL  is  achieved  in  this  pro¬ 
cess  by  making  sure  the  porosity  in  each  layer  approximately  equal 
to  the  prescribed  value.  Thus,  although  the  final  porosity  is  not 
exactly  equal  to  the  prescribed  value,  but  is  an  acceptable  approx¬ 
imation.  Fig.  2(b)  shows  the  microstructure  of  the  reconstructed 
carbon  paper  GDL  with  the  size  of  150  x  150  x  127  lattice  grid  and 
a  prescribed  porosity  of  0.791,  where  x  or  y  is  the  in-plane  direction 
and  z  the  through-plane  direction.  This  sample  is  also  considered  as 
the  original  sample  without  compression  and  PTFE  in  the  following 
study. 

The  compressed  samples  are  reconstructed  under  the  assump¬ 
tion  that  the  compression  process  only  change  the  thickness  of  the 


3.  Model  validation 


To  validate  the  present  model,  two  examples  of  single-phase 
flow  will  be  simulated:  (i)  fully-developed  flow  in  a  rectangular 
channel  and  (ii)  flow  through  body-centered  arrays  of  spheres.  For 
the  first  problem,  it  is  noted  that  the  analytical  solution  of  the  veloc¬ 
ity  profile  in  a  fully-developed  flow  in  a  channel  with  a  rectangular 
cross  section  having  width  2a  and  height  2b  is  [28] 


V(x,y)  =  V0 


fy\2  y^cosh(flx/fo)  cos(fey/b) 
\b)  2-^  cosh^a/h) 

i=i 


(21) 


with  V0  =  b2Apl(2pL)  and  =  (2/  —  1  )tt/2,  i  =  1,2, . . .,  n,  where  A p  is 
the  pressure  drop  over  the  channel  length  L,  and  p  is  the  dynamic 
viscosity  of  the  fluid.  Fig.  3(a)  shows  the  comparison  of  velocity 
profile  along  center  line  of  the  longer  side  obtained  from  the  ana¬ 
lytical  solution  given  by  Eq.  (21 )  with  prescribed  pressure  gradient 
A p/L  =  7.5  x  10-3  and  p  =  17.5,  and  the  LBM  simulation  for  two  rect¬ 
angular  channels  with  the  cross  section  size  of  60  x  30  and  80  x  30, 
respectively.  It  is  clearly  shown  that  the  simulation  results  are  in 
excellent  agreement  with  the  analytical  solution.  The  correspond¬ 
ing  pressure  drop  along  the  channel  is  presented  in  Fig.  3(b)  which 
shows  that  the  average  pressure  drops  linearly  along  the  channel 
length. 

The  validation  of  the  permeability  for  flow  through  body- 
centered  cubic  arrays  of  spheres  will  be  considered  next.  In  LBM 
simulation,  the  permeability  k  of  the  porous  medium  can  be  cal¬ 
culated  according  to  Darcy’s  law  under  the  low  Reynolds  number 
condition 

ud  =  J(Vp  +  F)  (22) 

where  ud  is  the  Darcy  velocity  as  the  superficial  velocity  of  the  fluid 
and  F  is  the  body  force  acting  on  the  fluid.  The  permeability  k  is 
generally  described  in  terms  of  the  porosity  0.  The  Kozeny-Carman 
(KC)  equation  is  the  most  widely  used  semi-empirical  relationship, 
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Fig.  3.  Comparison  of  the  LBM  simulation  results  with  the  analytical  solution  (solid 
line)  for  velocity  profile  of  fully-developed  flow  along  the  centerline  of  the  longer 
side  in  two  rectangular  channels,  (a)  Flow  profile  along  the  centerline  of  the  longer 
side  in  rectangular  channels  (with  cross  sections:  (H  )  60  x  30;  (•  )  80  x  30).  (b) 
Pressure  drop  along  the  channel. 


given  by 


k  =  C 


03 

(1-0)2 


(23) 


where  C  is  the  Kozeny-Carman  constant  and  C=d2/ 180  has  been 
used  for  packed-spheres  porous  media  with  sphere  diameter  d. 

Fig.  4  shows  the  comparison  of  simulation  results  of  the  dimen¬ 
sionless  permeability  k/d2  in  the  packed  sphere  porous  medium 
with  the  analytical  solutions  given  by  Zick  and  Homsy  [29]  and 
the  KC  relation  given  by  Eq.  (23).  It  can  be  seen  that  the  LBM 
simulation  results  are  in  good  agreement  with  the  analytical  solu¬ 
tion  for  the  whole  range  of  0,  but  the  KC  relation  shows  large 
error  as  the  porosity  increases  because  the  KC  equation  was  deter¬ 
mined  empirically  for  densely  packed  spheres  with  small  value 
of  0. 

After  this  validation,  the  viscosity  effect  on  the  permeability  will 
now  be  compared  based  on  the  MRT  and  BGK  models  for  the  recon¬ 
structed  carbon  paper  GDL  (as  shown  in  Fig.  2(b))  with  the  average 
porosity  of  0.791.  Fig.  5  shows  that  the  dimensionless  permeabil¬ 
ities  /</D2  calculated  by  the  MRT  model  are  constant  for  various 
viscosities,  while  the  permeability  calculated  in  BGK  model  is  sig¬ 
nificantly  increased  with  the  increase  of  viscosity.  Therefore,  the 
results  from  the  present  MRT  model  do  not  suffer  from  the  problem 


Fig.  4.  Comparison  of  the  LBM  simulated  permeability  with  the  analytical  solution 
and  the  KC  relation  for  body-centered  cubic  spheres  arrays. 


of  viscosity  dependent  permeability,  and  is  capable  of  simulating 
single-phase  flow  through  porous  media  more  accurately. 

4.  Results  and  discussion 

In  the  LBM  model,  the  simulation  parameters  are  in  terms  of  the 
lattice  units  that  do  not  have  physical  units.  To  relate  the  physical 
space  with  lattice  space,  a  length  scale  /0,  a  time  scale  t0  and  a  mass 
scale  m0  should  be  chosen  such  as  l0  =  1.5  x  10-6  m,  t0  =  5.0  x  10-7  s, 
and  m0  =  5.0  x  10-17  kg  in  present  study.  As  a  result,  the  physical 
parameters  of  velocity  up,  pressure  pp,  and  permeability  kp  can  be 
calculated  from  the  quantities  in  lattice  system  (subscripted  by  L) 
and  the  scale  parameters  as  follows: 

Up  =  Uilf,  Pp=PlT^’  kP  =  kLl0  (24) 

h)  to  fp 

For  the  reconstructed  carbon  paper  GDL,  the  fiber  diameter 
of  7.5  |jim  is  prescribed.  In  order  to  reduce  the  influence  of  the 
fiber  on  the  inlet/outlet  boundary,  five  more  lattices  are  added 
between  the  GDL  and  inlet/outlet  boundary.  For  in-plane  flow, 
pressure  boundaries  are  specified  on  the  inlet  and  outlet,  respec¬ 
tively  (perpendicular  faces  to  x-direction  in  Fig.  2(b)),  and  two 


0.0  0.1  0.2  0.3  0.4  0.5 

Viscosity  (LB  Unit) 

Fig.  5.  Comparison  of  the  permeability  computed  based  on  MRT/LBM  and  BGK/LBM 
models  with  various  viscosities. 
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Fig.  6.  Darcy  velocity  and  permeability  as  a  function  of  pressure  gradient  in  uncom¬ 
pressed  and  compressed  carbon  papers,  (a)  Darcy  velocity  vs.  pressure  gradient,  (b) 
Permeability  at  different  pressure  gradient. 


solid  boundaries  are  specified  on  the  top  and  bottom  surfaces  (per¬ 
pendicular  faces  to  z-direction  in  Fig.  2(b)),  keeping  the  others 
periodic  boundary  conditions.  For  through-plane  flow,  except  for 
the  pressure  boundaries  condition  on  the  inlet  and  outlet  (perpen¬ 
dicular  faces  to  z-direction  in  Fig.  2(b)),  the  rest  four  faces  are  all 
specified  to  be  periodic  boundary.  Additionally,  the  reproducibility 
of  the  reconstructed  GDL  must  be  ensured.  In  the  present  study, 
with  the  GDL  size  of  150  x  150  x  127  lattice  and  the  length  scale 
of  1.5  |jim,  the  GDLs  are  firstly  reconstructed  four  times  with  the 
same  porosity.  The  differences  among  the  calculated  permeabil¬ 
ity  of  the  four  GDLs  are  less  than  5%,  so  the  reconstructed  GDL 
size  (150  x  150  x  127  lattice)  can  meet  the  requirement  of  repro¬ 
ducibility,  and  the  corresponding  GDL  thickness  is  190.5  [xm  (i.e., 
127  x  1.5  |xm),  which  is  the  same  as  the  Tory  060  carbon  paper 
GDL. 

Fig.  6(a)  shows  the  relationship  between  the  Darcy  velocity  and 
the  pressure  gradient  for  the  original  and  25%  compressed  GDL 
based  on  MRT/LBM  simulation.  It  clearly  shows  that  the  Darcy 
velocity  is  linearly  proportional  to  the  pressure  gradient  for  both 
in-plane  and  through-plane  flows  within  the  calculated  pressure 
gradient  bound,  and  the  corresponding  permeability  is  almost  con¬ 
stant  with  pressure  gradient  as  shown  in  Fig.  6(b).  Therefore,  it  can 
be  said  that  the  Darcy  law  is  valid  and  the  inertial  force  can  be 
neglected  in  flow  through  carbon  paper  GDL  if  the  pressure  gradient 
does  not  exceed  the  bound  in  this  study. 


4.1  Velocity,  streamline  and  pressure  distribution 

Fig.  7  shows  the  3D  distribution  of  velocity,  streamline  and  pres¬ 
sure  in  the  carbon  paper  GDL  obtained  from  the  present  MRT/LBM 
model.  It  is  seen  that  the  flow  field  in  GDL  is  quite  complicated  for 
both  in-plane  and  through-plane  due  to  the  complicated  structure 
of  the  GDL.  In  each  slices,  the  magnitude  of  the  velocity  greatly 
depends  on  location  of  the  void  space.  It  can  be  seen  from  these 
graphs  that  the  main  flow  path  is  through  larger  pores  because  of 
their  small  flow  resistance,  and  it  can  be  expected  that  the  pressure 
would  drop  rapidly  at  the  locations  containing  dense  fibers. 

4.2.  Permeability 

It  is  observed  from  Fig.  7  that  the  streamlines  in  through-plane 
flow  are  a  little  more  tortuous  than  in-plane  flow.  This  implies  that 
liquid  particles  in  through-plane  flow  will  detour  more  obstacles 
resulting  in  larger  flow  resistance.  Thus,  the  in-plane  permeabil¬ 
ity  always  has  a  larger  value  than  through-plane  in  a  given  carbon 
paper  GDL  as  shown  in  Fig.  8. 


4.2.1.  Effects  of  compression 

The  effects  of  compression  ratio  (r)  on  in-plane  and  through- 
plane  permeabilities  of  carbon  papers  are  shown  in  Fig.  8,  where 
the  in-plane  permeability  of  Tory  060  obtained  by  Feser  et  al. 
[13]  and  the  through-plane  permeability  of  Tory  090  obtained  by 
Gostick  et  al.  [30]  are  also  plotted  for  comparison  purposes.  The 
in-plane  permeability  obtained  by  the  present  model  agrees  well 
with  measurements  by  Feser  et  al.  [13],  with  the  divergence  less 
than  15%  in  the  whole  range  of  compression  ratio.  On  the  other 
hand,  no  measurements  are  available  in  literature  for  the  through- 
plane  permeability  under  compression.  Flowever,  Gostick  et  al. 
[30]  measured  the  through-plane  permeability  for  uncompressed 
Toray  090  to  be  8.3  x  10-12  m2  while  Mathias  et  al.  [7]  measured 
the  through-plane  permeability  of  uncompressed  Toray  060  to  be 
approximately  in  the  range  of  5-10  x  10-12  m2.  The  simulated  value 
of  7.62541  x  10-12  m2  (with  no  compression  for  r  =  0)  falls  in  the 
middle  of  these  values. 

In  practice,  it  is  more  convenient  to  describe  the  permeability  as 
a  function  of  porosity,  such  as  KC  relation  given  by  Eq.  (23 ).  Since  the 
KC  relation  was  originally  determined  for  granular  porous  media 
and  is  applicable  only  for  low  porosity  materials  (as  shown  in  Fig.  4), 
therefore,  the  KC  relation  does  not  always  predict  the  correct  per¬ 
meability  for  fibrous  media  [31],  which  is  similar  to  carbon  paper 
GDL  with  high  porosities.  Flowever,  Feser  et  al.  [13]  found  that  KC 
relation  can  accurately  describe  the  permeability  of  carbon  paper 
in  a  narrow  porosity  range  by  fitting  a  proper  Kozeny-Carman  con¬ 
stant  C.  On  the  other  hand,  Tomadakis  and  Robertson  [32  ]  suggested 
that  the  permeability-porosity  KC  relationship  (Eq.  (23))  of  carbon 
paper  (which  has  randomly  overlapping  fiber  structure)  should  be 
replaced  by 


k  =  C 


4> 

(In  4>f 


(25) 


Tomadakis  and  Robertson  [32]  also  proposed  a  more  compre¬ 
hensive  relation  to  predict  the  anisotropic  permeability  of  in-plane 
and  through-plane,  based  on  the  A-base  method  of  Johnson  et  al. 
[33  ].  For  randomly  overlapping  fiber  structures,  the  A-base  relation 
is  described  as  [32] 


81n2</>(l  -<t>)a[{a  +  1)4>~4>p]2 


(26) 


where  R  is  the  radius  of  the  fiber,  0P  the  percolation  threshold  and 
a  a  constant  depending  on  the  structure  and  the  flow  direction. 
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Fig.  7.  Lattice  Boltzmann  simulation  results  for  single-phase  flow  through  3D  carbon  paper  GDL  (pressure  is  normalized  to  the  outlet  pressure),  (a)  Velocity  and  streamline 
for  through-plane  flow,  (b)  Pressure  contous  for  through-plane  flow,  (c)  Velocity  and  streamline  for  in-plane  flow,  (d)  Pressure  contous  for  in-plane  flow. 


Tomadakis  and  Robertson  [32]  found  that  for  in-plane  and  through- 
plane  flow  in  carbon  paper  where  the  percolation  threshold  </>p  is 
0.11,  the  constants  a  in  Eq.  (26)  are  0.521  and  0.785,  respectively. 
Gostick  et  al.  [14]  found  that  Eq.  (26)  with  a  =  0.521  fitted  the  in¬ 
plane  permeability  data  of  their  carbon  paper  with  fiber  diameter 
of  9.2  |jim  very  well. 

Fig.  9  shows  a  comparison  of  the  simulated  in-plane  and 
through-plane  permeabilities  of  carbon  paper  with  fiber  diame¬ 
ter  of  7.5  |jim  as  a  function  of  porosity  with  the  predicted  results 
from  the  relation  given  by  Eq.  (25)  and  A-base  relation  given  by 
Eq.  (26).  For  Eq.  (25),  the  fitted  constants  C  are  8.9504  x  10-13  and 
6.2805  x  10-13  for  in-plane  and  through-plane,  respectively.  From 
the  figure,  it  is  shown  that  Eq.  (25)  for  both  in-plane  and  through- 
plane  permeabilities  match  with  the  simulated  results  very  well. 
On  the  other  hand,  the  A-base  relation  is  shown  to  underestimate 
the  simulated  permeability  by  35%,  although  it  correctly  predicts 
the  trend  of  permeability  as  a  function  of  porosity. 

4.2.2.  Effects  ofPTFE  content 

Following  the  method  discussed  in  Section  2.3,  several  GDL 
samples  with  different  PFTE  contents  are  reconstructed.  The  simu¬ 
lated  in-plane  and  through-plane  permeabilities  are  presented  as  a 


Fig.  8.  Comparison  of  in-plane  and  through-plane  permeabilities  with  measure¬ 
ments:  (★ )  Feser  et  al.  [13];  (A )  Gostick  et  al.  [30]. 
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Porosity 


Fig.  9.  Comparison  of  the  LBM  simulated  in-plane  and  through-plane  permeabilities 
with  predictions  from  Eqs.  (25)  and  (26).  (a)  In-plane;  (b)  through-plane. 


function  of  porosity  in  Fig.  10.  However,  the  simulated  permeabil¬ 
ities  do  not  fit  Eq.  (25)  nor  the  A-base  relation  Eq.  (26)  very  well 
and  therefore  are  not  presented.  This  may  be  attributed  to  the  more 
complicated  structural  changes  after  the  PTFE  is  added.  It  is  obvious 


Porosity 


Fig.  10.  Comparison  of  the  simulation  permeability  with  predictions  from  Eq.  (27). 
Solid  line  for  in-plane  permeability  with  C=2.0  and  a  =  3.7681.  Dashed  line  for 
through-plane  permeability  with  C= 2.57  and  a  =  3.1355. 


from  Fig.  2(c)  that  the  carbon  paper  is  not  a  pure  fibrous  structure 
with  PTFE  added  due  to  the  formation  of  PTFE  films  and  blocks 
between  the  fibers.  Therefore,  a  more  general  empirical  formula 
[34]  is  used  to  fit  the  simulated  results: 

k  =  C0a  (27) 

where  C  and  a  are  two  constants.  It  is  shown  from  Fig.  10  that  Eq. 
(27)  with  the  fitting  values  of  C=  2.0  and  a  =  3.7681  for  in-plane  and 
C=  2.57  and  a  =  3.1355  for  through-plane,  matches  with  simulated 
values  well. 


4.3.  Tortuosity 


Generally,  tortuosity  is  defined  as  the  ratio  of  the  flow  path 
length  to  the  thickness  of  the  porous  media  along  the  flow  direc¬ 
tion.  It  is  used  to  calculate  the  effective  diffusivity  controlling  the 
mass  transfer  ability  of  reactant  in  GDL.  The  Bruggeman  equation 
[14]  has  been  commonly  used  for  tortuosity  prediction  in  PEMFC 
models,  which  is  given  by 


However,  Bruggeman  equation  was  determined  empirically 
from  isotropic  porous  media,  and  therefore  does  not  reflect  the 
anisotropic  property  of  carbon  paper  GDL.  Koponen  et  al.  [35] 
proposed  the  following  formula  for  in-plane  and  through-plane 
tortuosity  as 


X  —  1  Q 


1-0 

(0-0pf 


(29) 


where  is  the  percolation  threshold  while  a  and  m  are  constants. 

We  now  compute  the  tortuosity  of  the  carbon  paper  GDL  after 
the  flow  field  has  been  determined  using  the  MRT/LBM  model.  For 
this  purpose,  two  alternatives  can  be  used  to  obtain  the  tortuos¬ 
ity  by  computing  flow  path  lengths  numerically.  One  may  average 
over  the  actual  lengths  of  the  flow  lines  themselves,  and  the  other 
may  average  the  flow  lines  weighted  by  flux.  The  first  alternative  is 
suitable  for  molecular  diffusion  while  the  latter  alternative  is  more 
natural  for  fluid  flow  in  porous  media.  Both  of  these  tortuosities  of 
GDL  can  be  calculated  based  on  the  following  definition  formulas 
[36]: 


T,=  I 

(30a) 

XZjm/LMr,) 

T2  ££,«(*) 

(30b) 

where  N=  1000  is  the  number  of  the  streamlines  used  for  average; 
/(r2)  is  the  length  of  the  streamline  starting  at  r2-;  L  is  the  thickness 
of  the  porous  media  along  the  flow  direction;  u(  r2)  is  the  tangential 
velocity  at  the  starting  point  r,-;  t\  is  the  direct  average  tortuosity; 
t2  is  the  flux-weighted  tortuosity.  The  simulated  values  of  in-plane 
and  through-plane  tortuosities  are  shown  in  Fig.  11  as  a  function 
of  porosity.  It  can  be  seen  that  (i)  the  value  of  tortuosity  X\  is 
almost  the  same  as  tortuosity  r2,  although  X\  is  slightly  larger  espe¬ 
cially  at  low  porosity.  This  may  be  due  to  the  narrow  range  of  the 
streamline  length  at  high  porosity  of  GDL.  (ii)  The  through-plane 
tortuosity  is  higher  than  in-plane  (also  seen  in  Fig.  7),  and  (iii)  the 
Bruggeman  equation  (represented  by  a  solid  line)  falls  between  the 
simulated  in-plane  and  through-plane  tortuosities.  Matching  Eq. 
(29)  with  the  simulated  tortuosities  gives  a  =  0.326  and  m  =  0.863 
for  in-plane  flow  and  a  =  0.714  and  m  =  0.543  for  through-plane  flow, 
which  was  obtained  based  on  the  percolation  threshold  of  0.11  for 
carbon  paper  [32]. 
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Fig.  11.  The  calculated  tortuosity  (Eq.  (29))  as  a  function  of  porosity.  Short  dashed 
line  is  fitted  to  the  calculated  through-plane  points  by  Eq.  (30)  with  a  =  0.714  and 
m  =  0.543.  Dashed  dot  line  is  fitted  to  the  calculated  in-plane  points  by  Eq.  (30)  with 
a  =  0.326  and  m  =  0.863.  Solid  line  represents  the  Bruggeman  equation  (Eq.  (28)). 


4.4.  Fractal  model  for  permeability  prediction 


The  fractal  theory,  as  a  useful  tool  to  characterize  irregular  or 
disordered  objects,  has  been  applied  to  predict  the  permeability  of 
porous  media  and  GDL  [37-39].  In  these  studies,  the  through-plane 
permeability  is  expressed  as  a  function  of  area  fractal  dimension  Df 
and  tortuosity  fractal  dimension  Dt,  which  are  rather  difficult  to 
be  determined  using  the  box-counting  method  by  analyzing  SEM 
images  of  the  GDL. 

In  this  section,  the  fractal  model  will  be  used  to  predict  the  per¬ 
meability  of  carbon  paper  GDL  using  the  relationship  of  tortuosity 
and  porosity  obtained  in  Section  4.3.  Yu  and  Li  [38]  have  proposed 
a  relationship  between  porosity  and  area  fractal  dimension  Df  as 


<P  = 


dE—Df 


(31) 


where  df  is  Euclid  dimension;  Amin  and  Amax  are  the  minimum  and 
maximum  pore  sizes  in  GDL  with  their  ratio  assumed  to  be  10-2  in 
the  following  study. 

The  permeability  of  GDL  based  on  fractal  model  [37,39]  can  be 
expressed  as  follows: 


7rLl~DtX  Df 

12  8A  3  +  Dt-Df 


(32) 


where  L0  is  the  representative  length  along  the  flow  direction  and 
A  is  the  section  area  perpendicular  to  flow  direction  which  can  be 
written  as  [39] 


^max 


A  = 


1  —  (^min/^max) 


2  -Df 


Df 


40 


2 -Df 


(33) 


Substituting  Eqs.  (33)  and  (31 )  with  dE  =  2  into  Eq.  (32)  gives 


k  = 


T  1  —Dt  3  1  +Dt 
^0  Amax 


32 


2  -  Df  I 
3  +  Dt-Df  _ 


4> 

!-</> 


(34) 


On  the  other  hand,  the  flux-weighted  tortuosity  in  fractal  model 
can  be  obtained  as  follows: 


f-(L(X)/L0MX)mdX  3+Dt-Df  /  L0  \  Dt_1 
rmax  q(A)/(X)dA,  4  —  Df  Umax) 

JAmin 


(35) 


where /(A)  =  DfX^[nX~^+D  is  the  pore  number  probability  den¬ 
sity  with  diameter  A;  q( A)  =  (7rAp/128/xL(A))A4  is  the  flow  rate 
through  a  capillary  tube  of  diameter  A;  L(A)  =  L0A1_Dt  bases  on  the 
fractal  theory. 

Combining  Eqs.  (34)  and  (35),  the  tortuosity  fractal  dimension 
Dt  disappears,  thus  the  permeability  can  be  expressed  in  terms  of 
the  tortuosity  as 


X2max  2  -  Df  0 
32r  4  -  Df  1  -  0 


(36) 


where  the  maximum  pore  size  A,max  can  be  estimated  for  the  regular 
fiber  structure  for  through-plane  flow  [37]: 


Amax=  / - — - jd  ( 37 ) 

Y  ?t(i  -  \/4>) 

with  d  being  the  diameter  of  carbon  fiber.  Thus,  if  Eq.  (31)  is  used 
to  determine  Df  the  permeability  of  GDL  is  a  function  depending 
only  on  the  tortuosity  and  porosity,  which  is  deduced  from  the 
fractal  model.  The  predicted  through-plane  permeabilities  given 
by  Eq.  (36)  with  the  tortuosity  given  by  Eq.  (29)  with  a  =  0.714 
and  m  =  0.543  for  through-plane  is  plotted  in  Fig.  9(b).  As  shown 
from  this  figure,  the  through-plane  permeability  predicted  from 
the  fractal  model  is  very  close  to  KC  relation.  However,  the  frac¬ 
tal  model  described  above  cannot  accurately  predict  the  in-plane 
permeability  of  carbon  paper  GDL.  This  may  be  attributed  to  the 
non-representative  fractal  characteristics  of  the  cross  section  per¬ 
pendicular  to  the  in-plane  flow  direction. 


5.  Concluding  remarks 

In  this  work,  the  single-phase  flow  through  the  carbon 
paper  GDL  with  anisotropic  permeabilities  is  simulated  using 
the  MRT/LBM  model.  The  GDL  samples  are  reconstructed  using 
the  stochastic  method  by  considering  different  porosities  and 
microstructures  to  imitate  different  GDL  compression  ratios  and 
PTFE  contents.  In-plane  and  through-plane  permeabilities  as  well 
as  the  tortuosity  are  calculated  from  the  flow  field  by  simula¬ 
tion.  The  simulated  permeabilities  match  well  with  the  existing 
measurements  in  literature.  The  calculated  permeabilities  are  also 
compared  with  empirical  relationships,  which  predict  the  perme¬ 
ability  as  a  function  of  porosity  and  some  fitting  parameters  are 
determined.  In  addition,  the  tortuosity  of  the  carbon  paper  GDL 
obtained  from  simulation  also  shows  the  anisotropic  characteristic 
with  through-plane  tortuosity  higher  than  in-plane.  The  tortuosity 
is  used  in  a  fractal  model  to  predict  the  through-plane  permeabil¬ 
ity,  and  the  results  indicate  that  the  through-plane  permeability 
of  carbon  paper  GDL  based  on  the  fractal  model  and  the  KC  rela¬ 
tion  are  in  good  agreement  with  each  other.  Compared  with  the 
SRT/LBM  approach,  the  results  of  this  paper  show  that  the  MRT/LBM 
approach  has  the  advantages  of  numerically  more  stable  and  vis¬ 
cosity  independent  when  solid  boundary  is  present.  Therefore,  the 
MRT/LBM  approach  is  more  suitable  for  simulating  flow  in  porous 
media  at  the  pore  level.  Furthermore,  because  of  its  kinetic  nature, 
the  LBM  method  is  a  particularly  suitable  method  to  study  multi¬ 
phase  and  multicomponent  flow  problems  in  porous  media,  which 
are  currently  encountered  in  fuel  cell  researches.  Actually,  a  few 
complete  LBM  simulations  [40,41  ]  for  single-phase  multicompo¬ 
nent  flow  in  solid  oxide  fuel  cells  had  already  been  performed. 
To  our  knowledge,  however,  due  to  the  difficulties  of  simulating 
multicomponent  transport  coupling  with  two-phase  transport  in 
PEMFCs,  such  simulations  of  PEMFCs  are  not  presented  in  the  liter¬ 
ature  and  further  efforts  are  required  in  future. 
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